</div>
Tarea 3
Física Computacional
</div>
Autor: Celeste Castro Granados $\mathbb C \hat{e}l \mathbb s$
Primero llamaremos a las bibliotecas que vamos a necesitar:
using Plots
using LaTeXStrings
1. En clase vimos el método de Euler mejorado en el que nosotros introducimos el Jacobiano de forma analítica para ser evaluado por el método de integración. Realice una función que implemente el mismo método de integración, pero en lugar de calcular el Jacobiano con una función que lo evalúa analíticamente, ahora sea calculado de manera numérica utilizando la derivación numérica.¶
En general, una ecuación diferencial se representa por:
$$ \frac{d\vec x}{dt} = \vec g (\vec x, t) $$Y sabemos por lo visto en clase que el método de Euler mejorado realiza una aproximación mediante la derivada del sistema de ecuaciones utilizando la matriz jacobiana y el desarrollo de Taylor de la $\vec x (t)$ por medio de la implementación del siguiente algoritmo:
$$ \vec x_{n+1} = \vec x_n + \delta \vec g ( \vec x_n , t_n ) + \frac{1}{2} \delta^2 \mathbb{J}^* ( \vec x_n, t_n) \cdot \vec g^* ( \vec x_n, t_n ) $$En donde: $\mathbb{J}^* = \big ( \mathbb{J}, \frac{\partial \vec g}{\partial t} \big )$ con $\mathbb{J}$ el Jacobiano, y $\vec g^* = (\vec g, 1)$. (se agrega un 1 como entrada extra en $g^*$ para que no haya problemas al momento de multiplicar por $\mathbb{J}^*$, i.e. que las dimensiones concuerden)
Ahora, queremos una función que calcule de manera numérica la matriz $\mathbb{J}^*$. Para esto, recordamos que:
$$\mathbb{J}^* = \begin{bmatrix} \frac{\partial g_1}{\partial x_1} & \frac{\partial g_1}{\partial x_2} & ... & \frac{\partial g_1}{\partial x_n} & \frac{\partial g_1}{\partial t} \\ \frac{\partial g_2}{\partial x_1} & \frac{\partial g_2}{\partial x_2} & ... & \frac{\partial g_2}{\partial x_n} & \frac{\partial g_2}{\partial t} \\ . & & & \\ . & & & \\ . & & & \\ \frac{\partial g_n}{\partial x_1} & \frac{\partial g_n}{\partial x_2} & ... & \frac{\partial g_n}{\partial x_n} & \frac{\partial g_n}{\partial t} \end{bmatrix}$$Por tanto, necesitamos calcular las derivadas parciales de la función $\vec g$. Dichas parciales las calcularemos mediante la siguiente expresión:
$$ \partial_{x_j} g_{i,n} = \bigg (\frac{1}{2h} \bigg ) \bigg [ g_i \big (( x_{1,n}, ..., x_{n,n}, t_n ) + (0,...,h,...,0) \big ) - g_i \big (( x_{1,n}, ..., x_{n,n}, t_n ) - (0,...,h,...,0) \big ) \bigg] $$En donde, $h$ se encuentra en la j-ésima posición ya que estamos calculando la parcial respecto a $x_j$.
A continuación, procedemos a programar la función que nos de esta matriz:
function jacobiano_t(xi,h,f,t)
#xi-consiste en las variables del problema, donde calcularemos las parciales
#h-paso utilizado para obtener las derivadas parciales
#f-función a la cual sacaremos el jacobiano
#t-tiempo
#Definimos la matriz que vamos a ir llenando:
jacob=zeros(length(f(xi,t)),length(xi)+1)
for i in 1:length(f(xi,t)), j in 1:length(xi) #ciclo para ir llenando las entradas de la matriz
#xp1,xp2-vectores que tienen un aumento o disminución en la i-ésima entrda.
xp1=zeros(length(xi))
xp2=zeros(length(xi))
for k in 1:length(xi) #ciclo con el cual definiremos las entradas de los vectores donde evaluaremos la función para obtener la derivada parcial
if k==j #incremento a la j-ésima entrada del vector
xp1[k]=xi[j]+h
xp2[k]=xi[j]-h
else
xp1[k]=xi[j]
xp2[k]=xi[j]
end
end
jacob[i,j]=(f(xp1,t)[i]-f(xp2,t)[i])/(xp1[j]-xp2[j]) #entrada ij de la matriz, sin tomar en cuenta la entrada de la parcial temporal
jacob[i,length(xi)+1]=(f(xi,t+h)[i]-f(xi,t-h)[i])/(2h) #componentes de la columna de la parcial temporal
end
return jacob
end
Luego, como ya definimos la función que calcula de manera numérica la matriz $\mathbb{J}^*$, modificamos el algoritmo de Euler mejorado que hace uso del jacobiano con dependencia temporal visto en clase para integrarle nuestra función:
function Euler_mejorado_Jacobnum(edo,edo1,p_ini,t)
#edo-función de la ecuacion diferencial
#edo1-función de la ecuacion diferencial con un elemento 1 en la última entrada
#p_ini-condiciones iniciales de la posición
#t-tiempo
# J y g con dependencia temporal.
sol = zeros(length(t),length(p_ini))
sol[1,:] = p_ini
δ = t[2]-t[1]
for i in 1:(length(t)-1)
eval_edo = edo(sol[i,:],t[i])
g=edo1(sol[i,:],t[i]) #función por la que multiplicaremos el jacobiano
if length(eval_edo) == length(p_ini)
sol[i+1,:] .= sol[i,:] .+ δ .*eval_edo
else
sol[i+1,:] .= sol[i,:] .+ δ .*eval_edo[1:(end-1)]
end
sol[i+1,:] .+= 0.5*(δ^2) .*(jacobiano_t(sol[i,:],δ,edo,t[i]) *g )
end
return sol
end
2. Utilizando los integradores de Euler mejorado con Jacobiano analítico, Euler mejorado con Jacobiano numérico, Runge-Kutta de $2^°$ orden y Runge-Kutta de $4^°$ orden, encuentre un paso de integración que proporcione una precisión de $10^{-4}$ en la solución para la siguiente ecuación:¶
$$\dddot{y} + \ddot{y}^2 - 3\dot{y}^3 + cos^2 y = e^{-t} sin (3t)$$¶
Con condiciones iniciales $\ddot{y}(1)=1$, $\dot{y}(1)=2$, $y(1)=1$. En particular que sucede en el intervalo $t \in [1,2]$ y $t \in [1,2.1]$.¶
Recordamos primero los algoritmos vistos en clase de los métodos que utilizaremos:
Euler mejorado con Jacobiano analítico¶
function Euler_mejorado_J_t_p(edo,p_ini,t,jacob)
# J y g con dependencia temporal
sol = zeros(length(t),length(p_ini))
sol[1,:] = p_ini
δ = t[2]-t[1]
for i in 1:(length(t)-1)
eval_edo = edo(sol[i,:],t[i])
if length(eval_edo) == length(p_ini)
sol[i+1,:] .= sol[i,:] .+ δ .*eval_edo
else
sol[i+1,:] .= sol[i,:] .+ δ .*eval_edo[1:(end-1)]
end
sol[i+1,:] .+= 0.5*(δ^2) .*(jacob(sol[i,:],t[i]) *eval_edo )
end
return sol
end
Runge-Kutta de $2^°$ orden¶
function RK_2(edo,x_ini,t)
sol = zeros( length(t) , length(x_ini) )
sol[1,:] .= x_ini
δ = t[2]-t[1]
for i in 1:length(t)-1
k1 = sol[i,:] .+ 0.5*δ .*edo(sol[i,:],t[i])
sol[i+1,:] .= sol[i,:] .+ δ*edo(k1,t[i]+0.5*δ)
end
return sol
end
Runge-Kutta de $4^°$ orden¶
function RK_4(edo,x_ini,t)
sol = zeros( length(t) , length(x_ini) )
sol[1,:] .= x_ini
δ = t[2]-t[1]
for i in 1:length(t)-1
k1 = edo(sol[i,:],t[i])
k2 = edo(sol[i,:] .+ 0.5*δ.*k1 , t[i] + 0.5*δ)
k3 = edo(sol[i,:] .+ 0.5*δ.*k2 , t[i] + 0.5*δ)
k4 = edo(sol[i,:] .+ δ.*k2 , t[i] + δ)
sol[i+1,:] .= sol[i,:] .+ (δ/6.0).*(k1 .+ 2.0.*k2 .+ 2.0.*k3 .+ k4)
end
return sol
end
Euler mejorado con Jacobiano numérico¶
En este caso, utilizaremos el algoritmo desarrollado en el problema 1.
Definición de EDO y Jacobiano analítico¶
Vamos a analizar ahora la ecuación que tenemos:
Para facilitar las cosas, podemos ver a $y$ como la posición:
$$ \Rightarrow \dot{y}= v \quad (\text{velocidad}) $$$$ \Rightarrow \ddot{y}=\dot{v}= a \quad (\text{aceleración}) $$Por lo que, la ecuación diferencial se transforma en:
$$ \dot{a} + a^2 -3v^3 + cos^2 y = e^{-t} sin(3t) $$Y tenemos las siguientes 3 EDO's de 1er orden:
$$ \dot{a} = e^{-t} sin(3t) - a^2 +3v^3 - cos^2 y $$$$ \dot{v} = a $$$$ \dot{y} = v $$Y utilizando las condiciones iniciales dadas sabemos que: $a(1) =1$ , $v(1)=2$, $y(1)=1$.
Por lo tanto, podemos definir las funciones g's como:
$$ g_1 = e^{-t} sin(3t) - a^2 +3v^3 - cos^2 y $$$$ g_2 = a $$$$ g_3 = v $$Luego, calculamos de manera analítica el Jacobiano:
$$\mathbb{J}^* = \begin{bmatrix} \frac{\partial g_1}{\partial a} & \frac{\partial g_1}{\partial v} & \frac{\partial g_1}{\partial y} & \frac{\partial g_1}{\partial t} \\ \frac{\partial g_2}{\partial a} & \frac{\partial g_2}{\partial v} & \frac{\partial g_2}{\partial y} & \frac{\partial g_2}{\partial t} \\ \frac{\partial g_3}{\partial a} & \frac{\partial g_3}{\partial v} & \frac{\partial g_3}{\partial y} & \frac{\partial g_3}{\partial t} \end{bmatrix}$$$$ \Rightarrow \mathbb{J}^* = \begin{bmatrix} -2a & 9v^2 & -2 sin(y) cos(y) & e^{-t} (2 cos (3t) -sin (3t)) \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \end{bmatrix}$$Por lo tanto, ya tenemos todo listo. A continuación, procedemos a programar la ecuación diferencial y el Jacobiano:
function edo_y_1(variables,t)
#variables=(a,v,y)
g1=(exp(-t)*sin(3*t))-(variables[1]^2)+3*(variables[2]^3)-((cos(variables[3]))^2)
g2=variables[1]
g3=variables[2]
return [g1,g2,g3,1] #Se agrega un 1 para que las dimensiones coincidan y se pueda multiplicar por el jacobiano
end
function jacob_edo_y(variables,t)
Jacob = zeros(3,4) #definimos una matriz de ceros de 3x4 para poder llenarla
Jacob[1,1] = -2*variables[1]
Jacob[1,2] = 9*variables[2]^2
Jacob[1,3] = -2*sin(variables[3])*cos(variables[3])
Jacob[1,4] = exp(-t)*(2*cos(3*t)-sin(3*t))
Jacob[2,1] = 1
Jacob[3,2] = 1
return Jacob
end
Y recordamos que también debemos definir la función de la ecuación diferencial sin el 1 extra en el arreglo que entrega como resultado:
function edo_y(variables,t)
#variables=(a,v,y)
g1=(exp(-t)*sin(3*t))-(variables[1]^2)+3*(variables[2]^3)-((cos(variables[3]))^2)
g2=variables[1]
g3=variables[2]
return [g1,g2,g3]
end
Finalmente, definimos las condiciones iniciales para poder utilizarlas al momento de calcular las soluciones con cada uno de los métodos:
cond_inic_y=[1,2,1]; #condiciones iniciales
A continuación, vamos a obtener las soluciones con cada uno de los métodos mencionados arriba, encontrando en cada caso el paso adecuado para el cual obtenemos una precisión de $10^{-4}$, i.e. que obtenemos al menos 4 cifras significativas:
Euler con Jacobiano analítico
Vamos a ver que sucede primero en el intervalo [1,1.8]:
t_euler_a1= collect(1:0.001:1.8)
solucion_euler_a1 = Euler_mejorado_J_t_p(edo_y_1,cond_inic_y,t_euler_a1,jacob_edo_y)
Notamos que en el intervalo [1,1.8], las soluciones se encuentran bien definidas, sin embargo, veremos que si nos movemos a los intervalos [1,1.9], [1,2] ó [1,2.1] las soluciones ya no se encontrarán bien definidas:
t_euler_a2= collect(1:0.001:2.1)
solucion_euler_a2= Euler_mejorado_J_t_p(edo_y_1,cond_inic_y,t_euler_a2,jacob_edo_y)
Euler con Jacobiano numérico
Nuevamente veamos primero que sucede en el intervalo [1,1.8]:
t_euler_n1 = collect(1:0.00001:1.8)
solucion_euler_n1= Euler_mejorado_Jacobnum(edo_y,edo_y_1,cond_inic_y,t_euler_n1)
Y ahora para el intervalo [1,2.1]:
t_euler_n2= collect(1:0.00001:2.1)
solucion_euler_n2= Euler_mejorado_Jacobnum(edo_y,edo_y_1,cond_inic_y,t_euler_n2)
Runge-Kutta de orden 2
Primero para el intervalo [1,1.8]:
tRK2_1 = collect(1:0.0001:1.8)
solucion_RK2_1 = RK_2(edo_y,cond_inic_y,tRK2_1)
Y ahora para el intervalo [1,2.1]:
tRK2_2 = collect(1:0.0001:2.1)
solucion_RK2_2= RK_2(edo_y,cond_inic_y,tRK2_2)
Runge-Kutta de orden 4
Primero para el intervalo [1,1.8]:
tRK4_1 = collect(1:0.0001:1.8)
solucion_RK4_1 = RK_4(edo_y,cond_inic_y,tRK4_1)
Y ahora para el intervalo [1,2.1]:
tRK4_2 = collect(1:0.0001:2.1)
solucion_RK4_2 = RK_4(edo_y,cond_inic_y,tRK4_2)
En todas las soluciones que obtuvimos, observamos que todos los métodos tienen problemas a partir de aproximadamente 1.8, ya que ahí, se presentan indeterminaciones.
Finalmente, vamos a graficar la solución con el método de Rugen-Kutta de orden 4 y con el método de Euler con Jacobiano numérico para poder visualizar su comportamiento cerca de la zona de indeterminación que detectamos arriba, así como también, para comprobar que el algoritmo desarrollado en el problema 1 funciona correctamente:
plot(title="Solución en el intervalo [1, 1.8]",xlabel="t")
plot!(t_euler_n1,solucion_euler_n1[:,1],label="g_1(t) con Euler")
scatter!(tRK4_1,solucion_RK4_1[:,1],m=:cross,label="g_1(t) con RK4")
plot!(t_euler_n1,solucion_euler_n1[:,2],label="g_2(t) con Euler")
scatter!(tRK4_1,solucion_RK4_1[:,2],m=:cross,label="g_2(t) con RK4")
plot!(t_euler_n1,solucion_euler_n1[:,3],label="g_3(t) con Euler")
scatter!(tRK4_1,solucion_RK4_1[:,3],m=:cross,label="g_3(t) con RK4")
Primero, notamos que las soluciones con ambos métodos son similares, por lo que, podemos concluir que el algoritmo del problema 1 trabaja adecuadamente. Por otra parte, observando la gráfica vemos que efectivamente tenemos una variable que diverge y por tanto, no podemos obtener una solución acertada en dicha zona con ninguno de los métodos. De hecho, para asegurarnos aún más, podemos tomar nuestro intervalo de tiempo como la unión de dos intervalos, en donde el segundo intervalo va a ir desde 1.7 hasta 2.1 pero con muchos más puntos que en el primer intervalo para tratar de visualizar mejor el comportamiento justo en la zona donde habíamos visto que tenemos problemas:
t_euler_n21 = collect(1.7:0.08:2.1) # t entre 1.7 y 2.1
time=vcat(t_euler_n1,t_euler_n21)
solucion_analitico_21 = Euler_mejorado_J_t_p(edo_y_1,cond_inic_y,time,jacob_edo_y)
#solucion_numerico_21=Euler_mejorado_Jacobnum(edo_y,edo_y_1,cond_inic_y,time)
plot(title="Solución en el intervalo [1, 2.1]",xlabel="t")
plot!(time,solucion_analitico_21[:,1],label="g_1(t)")
plot!(time,solucion_analitico_21[:,2],label="g_2(t)")
plot!(time,solucion_analitico_21[:,3],label="g_3(t)")
Efectivamente hay problemas a partir de 1.75 aprox, y de hecho ni siquiera logra graficarlo bien.
3. Modifique las funciones integradoras que programamos en clase (Euler mejorado, Runge-Kutta de $2^°$ orden y Runge-Kutta de $4^°$ orden) para que se puedan pasar argumentos opcionales al sistema de ecuaciones diferenciales de tal manera que se puedan modificar parámetros al momento de hacer la integración. Una vez hecho esto, obtenga la solución de la siguiente ecuación diferencial:¶
$$ \ddot{x} + \frac{1}{10} \dot{x} + 4 sin x = \frac{1}{2} sin (\omega t) $$¶
Con condiciones iniciales $\dot{x}(0) = 1$ y $x(0) = 0$ y para $\omega = {1, 1.1, 1.2, 1.3, ..., 2.9, 3}$. Realice la gráfica de algunas soluciones que representen el comportamiento de la solución en función de $\omega$ y explique su resultado. La solución debe tener una precisión de por lo menos $10^{−5}$.¶
Euler mejorado¶
function Euler_mejorado_J_t_p_omega(edo,p_ini,t,jacob,ω)
#anexamos omega, que va a ser el parámetro que vamos a poder modificar
# J y g con dependencia temporal
sol = zeros(length(t),length(p_ini))
sol[1,:] = p_ini
δ = t[2]-t[1]
for i in 1:(length(t)-1)
eval_edo = edo(sol[i,:],t[i],ω) #Añadimos ω porque ahora la edo depende de ella
if length(eval_edo) == length(p_ini)
sol[i+1,:] .= sol[i,:] .+ δ .*eval_edo
else
sol[i+1,:] .= sol[i,:] .+ δ .*eval_edo[1:(end-1)]
end
sol[i+1,:] .+= 0.5*(δ^2) .*(jacob(sol[i,:],t[i],ω) *eval_edo )
end
return sol
end
Runge-Kutta de $2^°$ orden¶
function RK_2_omega(edo,x_ini,t,ω)
#anexamos ω, que va a ser el parámetro que vamos a poder modificar
sol = zeros( length(t) , length(x_ini) )
sol[1,:] .= x_ini
δ = t[2]-t[1]
for i in 1:length(t)-1
k1 = sol[i,:] .+ 0.5*δ .*edo(sol[i,:],t[i],ω)
sol[i+1,:] .= sol[i,:] .+ δ*edo(k1,t[i]+0.5*δ,ω)
end
return sol
end
Runge-Kutta de $4^°$ orden¶
function RK_4_omega(edo,x_ini,t,ω)
#anexamos ω, que va a ser el parámetro que vamos a poder modificar
sol = zeros( length(t) , length(x_ini) )
sol[1,:] .= x_ini
δ = t[2]-t[1]
for i in 1:length(t)-1
k1 = edo(sol[i,:],t[i],ω)
k2 = edo(sol[i,:] .+ 0.5*δ.*k1 , t[i] + 0.5*δ,ω)
k3 = edo(sol[i,:] .+ 0.5*δ.*k2 , t[i] + 0.5*δ,ω)
k4 = edo(sol[i,:] .+ δ.*k2 , t[i] + δ,ω)
sol[i+1,:] .= sol[i,:] .+ (δ/6.0).*(k1 .+ 2.0.*k2 .+ 2.0.*k3 .+ k4)
end
return sol
end
EDO y Jacobiano¶
Vamos a analizar ahora la ecuación diferencial que tenemos:
Tenemos:
$$ \ddot{x} + \frac{1}{10} \dot{x} + 4 sin x = \frac{1}{2} sin (\omega t) $$Con condiciones iniciales $\dot{x}(0) = 1$ y $x(0) = 0$ y para $\omega = {1, 1.1, 1.2, 1.3, ..., 2.9, 3}$.
Siguiendo el análisi realizado en la ecuación del problema 2, podemos ver a $x$ como la posición:
$$ \Rightarrow \dot{x}= v \quad (\text{velocidad}) $$$$ \Rightarrow \ddot{x}=\dot{v} \quad (\text{aceleración}) $$Por lo que, la ecuación diferencial se transforma en:
$$ \dot{v} + \frac{1}{10} v + 4 sin (x) = \frac{1}{2} sin (\omega t) $$Y tenemos las siguientes 2 EDO's de 1er orden:
$$ \dot{v} = \frac{1}{2} sin (\omega t) - \frac{1}{10} v - 4 sin (x) $$$$ \dot{x} = v $$Y utilizando las condiciones iniciales dadas sabemos que: $v(0)=1$, $x(0)=0$
Por lo tanto, podemos definir las funciones g's como:
$$ g_1=\frac{1}{2} sin (\omega t) - \frac{1}{10} v - 4 sin (x) $$$$ g_2 = v $$A continuación, procedemos a programar la ecuación diferencial:
function edo_x(variables,t,ω)
#variables=(v,x)
g1=(0.5*sin(ω*t))-((1/10)*variables[1]) - (4*sin(variables[2]))
g2=variables[1]
return [g1,g2,1]
end
Definimos las condiciones iniciales y el intervalo de tiempo que ocuparemos:
cond_inic_x=[1,0]
tiempo=collect(0:0.1:10*π);
Y calculamos su Jacobiano de manera analítica:
$$\mathbb{J}^* = \begin{bmatrix} \frac{\partial g_1}{\partial v} & \frac{\partial g_1}{\partial x} & \frac{\partial g_1}{\partial t} \\ \frac{\partial g_2}{\partial v} & \frac{\partial g_2}{\partial x} & \frac{\partial g_2}{\partial t} \\ \end{bmatrix}$$$$ \Rightarrow \mathbb{J}^* = \begin{bmatrix} -\frac{1}{10} & -4 cos(x) & \frac{1}{2}\omega cos(\omega t) \\ 1 & 0 & 0 \end{bmatrix}$$A continuación, procedemos a programar este Jacobiano:
function jacob_edo_x(variables,t,ω)
Jacob = zeros(2,3) #definimos una matriz de ceros de 2x3 para poder llenarla
Jacob[1,1] = -1/10
Jacob[1,2] = -4*cos(variables[2])
Jacob[1,3] = 0.5*ω*cos(ω*t)
Jacob[2,1] = 1
return Jacob
end
Y por último, recordamos que para los métodos de Runge-Kutta debemos volver a definir la EDO pero sin agregarle el 1 extra al arreglo que regresa la función.
function edo_x_RK(variables,t,ω)
#variables=(v,x)
g1=(0.5*sin(ω*t))-((1/10)*variables[1]) - (4*sin(variables[2]))
g2=variables[1]
return [g1,g2]
end
Con esto, ya tenemos todo lo necesario para obtener la solución con los 3 métodos. Ahora, el valor de $\omega$ no influye en la obtención de las cifras significativas ya que se trata de un parámetro constante, por tanto, basta con ver que se obtienen las 5 cifras significativas para una ω en particular. Vamos a obtener las soluciones para $\omega=1.1$:
solucion_euler_ω = Euler_mejorado_J_t_p_omega(edo_x,cond_inic_x,tiempo,jacob_edo_x,1.1) #solucion con el metodo de Euler
solucion_RK2_ω = RK_2_omega(edo_x_RK,cond_inic_x,tiempo,1.1) #solucion con el metodo de Runge-Kutta de orden 2
solucion_RK4_ω = RK_4_omega(edo_x_RK,cond_inic_x,tiempo,1.1) #solucion con el metodo de Runge-Kutta de orden 4
Observamos que para los 3 métodos se cumple que se tienen cinco cifras significativas como minímo. Por lo tanto, procedemos ahora sí a realizar las gráficas para presentar las soluciones con distintas ω:
ωs=collect(1:0.1:3)
plot(title="Solución para v(t) con distintos valores de ω (Euler)",xlabel="t",ylabel="v(t)")
for ω_i in ωs
solucion = Euler_mejorado_J_t_p_omega(edo_x,cond_inic_x,tiempo,jacob_edo_x,ω_i)
plot!(tiempo,solucion[:,1],label ="ω=$(ω_i)", lw=1.0)
end
plot!()
plot(title="Solución para x(t) con distintos valores de ω (Euler)",xlabel="t",ylabel="x(t)")
for ω_i in ωs
solucion = Euler_mejorado_J_t_p_omega(edo_x,cond_inic_x,tiempo,jacob_edo_x,ω_i)
plot!(tiempo,solucion[:,2],label ="ω=$(ω_i)", lw=1.0)
end
plot!()
plot(title="Solución para v(t) con distintos valores de ω (Runge-Kutta 2)",xlabel="t",ylabel="v(t)")
for ω_i in ωs
solucion = RK_2_omega(edo_x_RK,cond_inic_x,tiempo,ω_i)
plot!(tiempo,solucion[:,1],label ="ω=$(ω_i)", lw=1.0)
end
plot!()
plot(title="Solución para x(t) con distintos valores de ω (Runge-Kutta2)",xlabel="t",ylabel="x(t)")
for ω_i in ωs
solucion = RK_2_omega(edo_x_RK,cond_inic_x,tiempo,ω_i)
plot!(tiempo,solucion[:,2],label ="ω=$(ω_i)", lw=1.0)
end
plot!()
plot(title="Solución para v(t) con distintos valores de ω (Runge-Kutta 4)",xlabel="t",ylabel="v(t)")
for ω_i in ωs
solucion = RK_4_omega(edo_x_RK,cond_inic_x,tiempo,ω_i)
plot!(tiempo,solucion[:,1],label ="ω=$(ω_i)", lw=1.0)
end
plot!()
plot(title="Solución para x(t) con distintos valores de ω (Runge-Kutta4)",xlabel="t",ylabel="x(t)")
for ω_i in ωs
solucion = RK_4_omega(edo_x_RK,cond_inic_x,tiempo,ω_i)
plot!(tiempo,solucion[:,2],label ="ω=$(ω_i)", lw=1.0)
end
plot!()
Podemos notar que tanto para $x(t)$, como $v(t)$, tenemos un fenómeno de resonancia. Asimismo, conforme el sistema va avanzando en el tiempo, la amplitud de las soluciones va creciendo, haciendo que mientras mayor sea el valor de $\omega$ también sea mayor el valor de la amplitud.
4. Considere el siguiente sistema de ecuaciones diferenciales:¶
$$ \frac{dx}{dt} = 10 (y-x) $$¶
$$ \frac{dy}{dt} = x (28-z)-y $$¶
$$ \frac{dz}{dt} = xy - \frac{8}{3} z $$¶
Realice una integración numérica del sistema de ecuaciones con el método de Runge-Kutta de $4^°$ y realice los siguientes ejercicios:¶
a) Elija una condición inicial aleatoria en el $x, y, z \in [−5, 5]$ y realice la gráfica de la solución en el plano xy, xz y yz y comente sus observaciones.¶
b) Ahora tome 3 soluciones parecidas. Para ello proponga 3 condiciones iniciales tales que en la coordenada en $y$ las condiciones iniciales difieran entre si en $10^{−2}$. Explique lo que observa.¶
Nota: De preferencia realice integraciones con tiempos relativamente largos en los que pueda apreciar diferencias en el comportamiento de las soluciones. Además asegure que cada una de las soluciones que encuentre sea convergente en por lo menos 5 cifras ($10^{−5}$), es decir encuentre un paso adecuado para realizar la integración.¶
Vamos a definir primero nuestro sistema de ecuaciones:
function edo_xyz(variables,t)
#variables(x,y,z)
dx=10*(variables[2]-variables[1])
dy=variables[1]*(28-variables[3])-variables[2]
dz=variables[1]*variables[2]-((8/3)*variables[3])
return[dx,dy,dz]
end
a) Utilizaremos la función rand para elegir una condición inicial aleatoria:¶
cond_inic_xyz=[-5rand(),5rand(),5rand()];
Gráficas de la solución:
plot(title="Solución de la EDO (Plano xy)",xlabel="x",ylabel="y")
tiempo = collect(0:0.01:60π)
sol_num_rk4 = RK_4(edo_xyz,cond_inic_xyz,tiempo)
plot!(sol_num_rk4[:,1],sol_num_rk4[:,2])
plot!(legend=false)
plot(title="Solución de la EDO (Plano xz)",xlabel="x",ylabel="z")
tiempo = collect(0:0.01:60π)
sol_num_rk4 = RK_4(edo_xyz,cond_inic_xyz,tiempo)
plot!(sol_num_rk4[:,1],sol_num_rk4[:,3])
plot!(legend=false)
plot(title="Solución de la EDO (Plano yz)",xlabel="y",ylabel="z")
tiempo = collect(0:0.01:60π)
sol_num_rk4 = RK_4(edo_xyz,cond_inic_xyz,tiempo)
plot!(sol_num_rk4[:,2],sol_num_rk4[:,3])
plot!(legend=false)
En las 3 gráficas obtenemos la misma forma pero vista desde distintos ángulos. Asimismo, podemos ver que conforme el tiempo va creciendo, aparecen más líneas que siguen la misma trayectoria.
Finalmente, vamos a comprobar que la solución tiene al menos 5 cifras significativas:
sol_num_rk4 = RK_4(edo_xyz,cond_inic_xyz,tiempo)
Efectivamente obtuvimos una solución con al menos 5 cifras significativas.
b) Ahora, para tomar 3 condiciones iniciales aleatorias y parecidas realizaremos lo siguiente:¶
#Partimos de la condicion inicial que ya teníamos:
cond_inic_1=[cond_inic_xyz[1],cond_inic_xyz[2],cond_inic_xyz[3]]
#Y para tener otras 2 condiciones iniciales cercanas, le vamos a sumar y restar 0.01 para que la diferencia entre ellas sea de 10^{-2}
cond_inic_2=[cond_inic_xyz[1],cond_inic_xyz[2]+0.01,cond_inic_xyz[3]]
cond_inic_3=[cond_inic_xyz[1],cond_inic_xyz[2]-0.01,cond_inic_xyz[3]]
plot(title="3 Soluciones parecidas de la EDO (Plano xy)",xlabel="x",ylabel="y")
tiempo = collect(0:0.01:60π)
solnum_rk4_1 = RK_4(edo_xyz,cond_inic_1,tiempo)
solnum_rk4_2 = RK_4(edo_xyz,cond_inic_2,tiempo)
solnum_rk4_3 = RK_4(edo_xyz,cond_inic_3,tiempo)
plot!(solnum_rk4_1[:,1],solnum_rk4_1[:,2], label="Cond. inicial 1")
plot!(solnum_rk4_2[:,1],solnum_rk4_2[:,2], label="Cond. inicial 2")
plot!(solnum_rk4_3[:,1],solnum_rk4_3[:,2], label="Cond. inicial 3")
plot!(legend=true, grid=true)
plot(title="3 Soluciones parecidas de la EDO (Plano xz)",xlabel="x",ylabel="y")
tiempo = collect(0:0.01:60π)
solnum_rk4_1 = RK_4(edo_xyz,cond_inic_1,tiempo)
solnum_rk4_2 = RK_4(edo_xyz,cond_inic_2,tiempo)
solnum_rk4_3 = RK_4(edo_xyz,cond_inic_3,tiempo)
plot!(solnum_rk4_1[:,1],solnum_rk4_1[:,3], label="Cond. inicial 1")
plot!(solnum_rk4_2[:,1],solnum_rk4_2[:,3], label="Cond. inicial 2")
plot!(solnum_rk4_3[:,1],solnum_rk4_3[:,3], label="Cond. inicial 3")
plot!(legend=true, grid=true)
plot(title="3 Soluciones parecidas de la EDO (Plano yz)",xlabel="x",ylabel="y")
tiempo = collect(0:0.01:60π)
solnum_rk4_1 = RK_4(edo_xyz,cond_inic_1,tiempo)
solnum_rk4_2 = RK_4(edo_xyz,cond_inic_2,tiempo)
solnum_rk4_3 = RK_4(edo_xyz,cond_inic_3,tiempo)
plot!(solnum_rk4_1[:,2],solnum_rk4_1[:,3], label="Cond. inicial 1")
plot!(solnum_rk4_2[:,2],solnum_rk4_2[:,3], label="Cond. inicial 2")
plot!(solnum_rk4_3[:,2],solnum_rk4_3[:,3], label="Cond. inicial 3")
plot!(legend=true, grid=true)
Observando las gráficas anteriores, podemos notar que todas las soluciones tienen la misma forma sin importar la condición inicial, sin embargo si se puede apreciar un ligero desfase entre ellas (ya que se pueden distinguir los 3 colores), por lo que no podemos decir que son iguales, lo cual tiene sentido ya que las soluciones dependen justo de las condiciones iniciales que se tengan.
Finalmente, vamos a comprobar que para cada una de las soluciones tengamos al menos 5 cifras significativas:
solnum_rk4_1 = RK_4(edo_xyz,cond_inic_1,tiempo)
solnum_rk4_2 = RK_4(edo_xyz,cond_inic_2,tiempo)
solnum_rk4_3 = RK_4(edo_xyz,cond_inic_3,tiempo)
Efectivamente, para los 3 casos tenemos una solución con al menos 5 cifras significativas.